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Abstract. In real- world geophysical applications (such as predicting the climate change), the 
reduced models of real-world complex multiscale dynamics are used to predict the response of the 
actual multiscale climate to changes in various global atmospheric and oceanic parameters. However, 
while a reduced model may be adjusted to match a particular dynamical regime of a multiscale pro- 
cess, it is unclear why it should respond to external perturbations in the same way as the underlying 
multiscale process itself. In the current work, the authors study the statistical behavior of a reduced 
model of the linearly coupled multiscale Lorenz 96 system in the vicinity of a chosen dynamical 
regime by perturbing the reduced model via a set of forcing parameters and observing the response 
of the reduced model to these external perturbations. Comparisons are made to the response of the 
underlying multiscale dynamics to the same set of perturbations. Additionally, practical viability 
of linear response approximation via the Fluctuation-Dissipation theorem is studied for the reduced 
model. 
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1. Introduction 

A reduced model for sfow variables of multiscale dynamics is a lower-dimensional 
dynamical system, which "resolves" (that is, qualitatively approximates in some ap- 
propriate sense) major large scale slow variables of the underlying higher-dimensional 
multiscale dynamics while at the same time being relatively simple and compu- 
tationally inexpensive to work with. This is important in real-world applications 
of contemporary science, such as geophysical science and climate change predic- 
tion |14II15II19II21II24II28II37| , where the actual underlying physical process is impossible 
to model directly, and its reduced approximation has to be designed for such a purpose. 
Reduced dynamics were used to model global circulation patterns |131ll8,"36''48[l51j. 
and large-scale features of tropical convection [23l[30] . Typically, reduced models of 
multiscale dynamics consist of simplified lower-dimensional dynamics of the original 
multiscale dynamics for the resolved variables, with additional terms and parame- 
ters which serve as replacements to the missing coupling terms with the unresolved 
variables of the underlying physical process. These extra parameters in the reduced 
model are usually computed to match a particular dynamical regime of the underlying 
multiscale dynamics [t|l^[71 [T^[TCl[T71[^ [nHM l Hg] . In particular, if the underlying 
multiscale process changes its dynamical regime (for example, in response to changes 
in its own forcing parameters), then the parameters of the corresponding reduced 
model have to be appropriately readjusted to match its dynamical regime to the new 
regime of the multiscale dynamics. 

In some real-world applications, such as the climate change prediction, the re- 
duced models of complex multiscale climate dynamics are used to predict the response 
of the actual multiscale climate to changes in various global atmospheric and oceanic 
parameters. However, while a reduced model may be manually adjusted to match a 
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particular dynamical regime of a multiscale process, it is unclear whether it should 
respond to identical external perturbations a priori in the same way as the multiscale 
process, without any extra readjustments. How do reduced models of multiscale dy- 
namics, adjusted to a particular dynamical regime, respond to external perturbations 
which force them out of this regime? Is their response similar to the response of the 
underlying multiscale dynamics to the same external perturbation? It is quite clear 
that the reduced dynamics evolve on a set with lower dimension than that of the full 
multiscale dynamics. How do the properties of this limiting set respond to changing 
external forcing parameter, in comparison to the full multiscale attractor? 

Here we develop a set of criteria for similarity of the response to small external per- 
turbations between slow variables of multiscale dynamics and those of a reduced model 
for slow variables only, determined through statistical properties of the unperturbed 
dynamics. We also carry out a computational study of the difference in responses of 
the full multiscale and deterministic reduced dynamics of the linearly coupled rescaled 
Lorenz '96 model from |4][7j to identical external perturbations. We compare and con- 
trast both the actual ( "ideal" ) responses of the multiscale and reduced models directly 
to finitely small perturbations of external forcing, and the linear response predictions 
of the reduced models via the Fluctuation-Dissipation theorem [TH3ll8l- [Tl1l27ll40] . Two 
different types of forcing perturbations are used: the time-independent Heaviside forc- 
ing, and the simple time-dependent ramp forcing. The manuscript is structured as 
follows. In Section[5]we formulate the standard averaging formalism to obtain the av- 
eraged slow dynamics from a general two-scale dynamical system. Section [3] describes 
statistically tractable criteria to ensure similarity of responses between a two-scale 
system and its averaged slow dynamics. In Section |4] we describe the first-order re- 
duced model approximation to a two-scale dynamics with linear coupling between the 
slow and fast variables, previously developed in [4l[6j[7]. In Section [5] we introduce 
the two-scale Lorenz '96 toy model which will be our testbed for this method. In 
Section [6] we present comparisons of the large features of the multiscale and reduced 
systems, including statistical comparisons as well as the ability of the reduced model 
to capture perturbation response of the multiscale system. Section [7] summarizes the 
results and suggests future work. 

2. Averaged slow dynamics for a general two-scale system 

A general two-scale dynamical system with slow variables x and fast variables y 
is usually represented as 



where, x{t) gM.^^ are the slow variables of the system, y{t)EM.^y are the fast vari- 
ables, and F and G are nonlinear differentiable functions. The integer parameters 
Nx <^ Ny are the dimensions of the slow and fast variable subspaces, respectively. 
Usually, a time-scale separation parameter is used to denote the difference in time 
scales between the slow and fast variables, however, here we omit it, as the framework 
for reduced models from [4l|6ll7| , which we use here, does not require such a parameter 
to be explicitly present. 

Under the assumption of "infinitely fast" y-variables, one can use the averaging 




(2.1) 
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formalism [3H1IM1I1S1S3 write the averaged system for sfow variables alone: 
dec - - f 

— F{x)=J F{x,y)d^i4y), (2.2) 

where fx^ is the invariant distribution measure of the fast limiting system 

^ = G(x,z), (2.3) 

with X above in (|2.3|) being a constant parameter. We express the slow solutions of the 
two-scale system in (12.11) and the averaged system in (12. 2p in terms of differentiable 
flows: 

x{t) = (^*(a;o,yo) fo'" ^^e two-scale system, (2.4a) 



XA{t) — (f)\{xQ) for the averaged system. (2-4b) 

It can be shown (see [3H1I3S1I1S1I1I1 and references therein) that if the time scale 
separation between x and y is large enough, then, for the identical initial conditions 
Xq and generic choice of yg i the solution xa {t) of the averaged system in (|2.2I) remains 
near the solution x{t) of the original two-scale system in (|2.ip for finitely long time. 

3. Criteria of similarity of responses to small external perturbations 
for general two-scale system and its averaged slow dynamics 

Let /i and fj,A denote the invariant distribution measures for the two-scale system 
in ()2.1|) and the averaged system in p.2|) . respectively. Also, let h{x) be a differen- 
tiable test function. Then, the statistically average values of h for both two-scale and 
averaged systems are given via 

(h) = / h{x)dn{x,y), (3.1a) 



{h)A^ Hx)d^lA{x). (3.1b) 



Now, consider the two-scale system in ()2.1|) . and the averaged system in (|2.2|) . per- 
turbed at the slow variables by a small time-dependent forcing Sf{t): 

^= F{x,y)+Sf{t), 

dy ^^-^^^ 



— ^F{x) + Sf{t). (3.2b) 

Then, for small enough df{t), the average responses 6{h){t) and 5{h)A{t) for the 
two-scale system in (|2.1I) and the averaged system in (|2.2I) . respectively, can be ap- 
proximated by the following linear response relations: 

6{h){t) = J^R{t-s)Sf{s)ds, R{t)=J^!^^f^^l^dfi{x,y), (3.3a) 
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d{h)Ait)^j\Ait~s)Sf{s)ds, R^{t)^J^^^^^^df,A{x). (3.3b) 

For details, see P^H51I MTTll44) . Above, it is clear that any differences between d{h){t) 
and 6{h)fj{t) are due to differences between R{t) and RA{t), since Sf is identical 
in both cases. The differences between R{t) and RA{t) are, in turn, caused by the 
differences between the flows </>* and (j^A, and the differences between the invariant 
distribution measures /i and ^a, which are difficult to quantify in practice. In what 
follows we express the differences between R{t) and RA{t) via statistically tractable 
quantities. First, we assume that the invariant measures /i and fj,A are absolutely 
continuous with respect to the Lebesgue measure, with distribution densities p{x,y) 
and pa{x), respectively: 

Rit) = J'-!^^^^^pix,y)dxdy, (3.4a) 

RAit)-J'-^^^PAix)dx. (3.4b) 

While it is known that purely deterministic processes may not have Lebesgue- 
continuity of their invariant measures [42l|43l[50], however, even small amounts of 
random noise, which is always present in real- world complex geophysical dynamics, 
usually ensure the existence of the distribution density. Integration by parts yields 

R{t)^- jh{<t,\x,y))^^^^dxdy, (3.5a) 

RA{t)^~ jh{^\{x))^P^dx. (3.5b) 

At this point, let us express p{x,y) as the product of its marginal distribution p{x)^ 
defined as 

P{x)^ j p{x,y)dy, (3.6) 
and conditional distribution p{y\x), given by 

^(^"")^7(^- ^'-'^ 

It is easy to check that the conditional distribution p{y\x) satisfies the identity 

j p{y\x)dy=l foraUa;. (3.8) 
Now, the formula for the linear response operator R{t) above can be written as 

R{t) = -j h{^*{x,y))p{y\x)^^dydx- j h{^\x,y))^^^p{x)dydx. (3.9) 
We now denote 



e*{x,y) = (f)\x,y)~(j)\{x), 



(3.10) 
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where £*{x,y) is small compared to either (f)*{x,y) or (p^ix) for relevant values of t, 
X and y. Then, for the second integral in the right-hand side of p.9|) we write 



(3.11) 



- j h{<t>\x^y)f-P^p{x)dyAx = - j (^|^^dy) W)p»d=r- 

- jvh{4>\{x))e\x,y)^P^-p{x)dydx = 0{\\e\\] 

where the first integral in the right-hand side is zero due to the condition in p.Sp . 
Neglecting the 0(||e||) terms in (|3.5p . we write 

R{t) = -j h{<j>'ix,y))p{y\x)^^dydx, (3.12a) 

RA{t)^~ j h{<p\{x))^^^dx. (3.12b) 

At this point, we express p{x) and pa{x) as exponentials 

p(a;) = e-^(^), pa{x) = e'^^^^l (3.13) 

where b{x) and bA{x) are smooth functions, growing to infinity as x becomes infinite. 
The latter yields 



R{t)^ J h{^*ix,y))^^pix,y)dydx, 



(3.14a) 



RA{t) = Jh{<j,\{x))^^^PA{x)dx. (3.14b) 

Replacing invariant measure averages with long-term time averages yields the follow- 
ing time correlation functions: 

R{t)=\im- h{x{s + t)) — {x{s))ds, (3.15a) 
RA{t)^ lim- hixAis + t))—^{xA{s))ds. (3.15b) 

r ->oo r Jq ox 

Taking into account the arbitrariness of h, we conclude that, in order for i?^(i) to 
approximate R{t) despite the fact that, for long times s, xa{s) diverges from x(s) even 
for identical initial conditions, we generally need three conditions to be approximately 
satisfied: 

1. For identical initial conditions, XA{t) should approximate x{t) (that is, 
e*{x,y) in (|3.10p should indeed be small) on the finite time scale of decay 
of the correlation functions in p.lS^ : 

2. bA{x) should approximate b{x), which means that the invariant distribution 
Pa{x) of the averaged system in (|2.2p should be similar to the a;-marginal 
p{x) of the invariant distribution of the two-scale dynamical system in (j2.ip : 
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3. The time autocorrelation functions of the averaged system in (|2.2[) should 
be similar to the time autocorrelation functions of the slow variables of the 
two-scale system in (|2.1[) . 
As a side note, observe that the nature of dependence of the conditional distribution 
p{y\x) on X does not play any role in the criteria for the similarity of responses. In 
particular, the exact factorization of p{x,y) into its x- and y-marginals (which means 
that p{y\x) is independent of a;) is not required, unlike what was suggested in 29 for 
the Gaussian invariant states. 

4. Practical implementation of the reduced model for a two-scale pro- 
cess with linear coupling 

As formulated above in Sections [5] and |31 the criteria of the response similarity 
are applicable for a broad range of dynamical systems with general forms of coupling 
and their averaged slow dynamics. However, the practical computation of the reduced 
model approximation to averaged slow dynamics depends on the form of coupling in 
the two-scale system [4j[6j[7]. In this work, we consider the linear coupling between 
the slow and fast variables in the two-scale system (|2.1|) . The linear coupling is the 
most basic form of coupling in physical processes, however, because of that it is also 
probably the most common form of coupling. For the linear coupling, the reduced 
model is constructed according to the method developed previously in ^ , which we 
briefly sketch below. 

We consider the special setting of (|2.1[) with linear coupling between x and y: 



dec 

— = f{x) + Lyy, 
^= g{y) + L^x, 



(4.1) 



where / and g are nonlinear differentiable functions, and and Ly are constant 
matrices of appropriate sizes. The corresponding averaged dynamics for slow variables 
from p.2p simplifies to 

f{x) + Lyzix), (4.2) 



dx 
'dt 

where z{x) is the statistical mean state of the fast limiting system 



^^g{z) + L,x, (4.3) 

with X treated as constant parameter. In general, the exact dependence of z{x) on x 
is unknown, except for a few special cases like the Ornstein-Uhlenbeck process [45] . 
Here, like in [4j[7], we approximate z{x) via the linear expansion 

z{x)^z* +CL^{x-x*), (4.4) 

where x* is the statistical average state of the full multiscale system in (14.11) . and 
z* = z{x*). The constant matrix C is computed as the time integral of the correlation 
function 



e^( f C(s)ds]c-\0), C(s)=\im-f z{t + s)z^(t)dt, 



(4.5) 



where z{t) is the solution of (|4.3p for x = x*. The above formula constitutes the 
quasi-Gaussian approximation to the linear response of z to small constant forcing 
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perturbations in (I4.3p . and is a good approximation when the dynamics in (j4.3p are 
strongly chaotic and rapidly mixing [IH?ll5 HTT]|?7[Hn] . With (|4.5p . the reduced system 
in (j4.2l) becomes the explicitly defined deterministic reduced model for slow variables 
alone: 

d T' 

— = /(a;) + L,r+i(a;-a;*), (4.6) 

where L^LyCL^- In what follows, the "zero-order" model refers to (|4.6|) with the 
last term set to zero (such that the coupling is parameterized only by the constant 
term LyZ*). For details, see [Ul^ and references therein. 

5. Testbed the rescaled Lorenz '96 system 

In the current work, we test the response of the reduced model for slow variables 
on the rescaled Lorenz '96 system with linear coupling [4], which is obtained from 
the original two-scale Lorenz '96 system |26| by appropriately rescaling dynamical 
variables to approximately set their mean states and variances to zero and one, re- 
spectively. Below we present a brief exposition of how the rescaled Lorenz '96 model 
is derived. 

5.1. The original two-scale Lorenz '96 system 

The original two-scale Lorenz '96 system [26, is given by 

A ■' 

±1 = Xi^i{xi+i - Xi-2) -x., + F^ ^ ^ 2/i J , 

j=i (5.1) 

^ m,j = - [V'i.j+i {yi,j-i - yi.j+2) - Vt.j +Fy+ x^Xi] , 

where 1 < i < N^, ^<j<J, and periodic boundary conditions Xi+ — Xi , yi+ j = yi,j 
and yi.j+j — Vi+ij- Here and Fy are constant forcing terms, A^, and Xy constant 
coupling parameters, and e is the time scale separation parameter. Throughout this 
paper we will consider systems with twenty slow variables (N^ = 20) and eighty fast 
variables (TV^ = 80, J = 4) . 

In Lorenz's original formulation [26 studying predictability in atmospheric-type 
systems, he begins with the uncoupled system 

x^^x^^i{xi+i-Xi^2)-Xi + F, i = l,...,N, (5.2) 

with periodic boundary conditions. This system has generic features of geophysical 
flows, namely a nonlinear advection-like term, linearly unstable waves, damping, forc- 
ing, mixing, and chaos |35| . The simple formulation, with invariance under index 
translation and a uniform forcing term F, allows for straightforward analysis - in 
particular the long-time statistics of each variable should be identical and will only 
depend on F. Additionally, the chaos and mixing of the system are simply regu- 
lated by the forcing, with decaying solutions for F near zero, periodic solutions for 
F slightly larger, weakly chaotic quasi-periodic solutions around ^^ = 6, and chaotic 
and strongly mixing systems around _F=16 and higher. Lorenz's two-time coupled 
system was introduced to study predictability and Lyapunov exponents of systems 
with subgrid phenomena on faster timescales, and one of the authors of the current 
work has recent results showing that coupling two chaotic systems can suppress chaos 
in the slower system [S]. 
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5.2. Rescaled Lorenz '96 model 

To simplify the analysis of coupling trends for the two-time system, we will scale 
out the dependence of the mean state and mean energy on the forcing term F. Due 
to the translational invariance, the long-term mean x and standard deviation a for 
the uncoupled system are the same for all Xi. So we rescale x and t as 

Xi — X I CTX'i ^ t — ^ 

a 

where the new variables x have zero mean and unit standard deviation, while their 
time autocorrelation functions have normalized scaling across different dynamical 
regimes (that is, different forcings F) for short correlation times. This rescaling 
was previously used in 12 7j . In the rescaled variables, the uncoupled Lorenz model 
becomes 



^ X \ , ^ ^ , X i F X / r' o N 

Xi_i + - \{Xi+i-Xi^2) \ (5.3) 



where x and a are functions of F. 

We similarly rescale the coupled two-scale Lorenz '96 model: 

dxi f , x\ Xi F^-x A.y>A 

I [Xi+i - Xi„2 ) 1 2 T / . 



at \ ax J Ox Ox J . , 

J=i (5.4) 

£—rr^ yt,3+i + — [yi,j-i-yz,j+2) — -+ . +^^^»^ 

dt \ ' Oy J ay Oy^ 



where {x, Ox} and {y, Oy} are the long term means and standard deviations of the 
uncoupled systems with Fx or Fy as constant forcing, respectively. It is this rescaled 
coupled Lorenz '96 system that we focus on for the closure approximation. 

Before any numerical tests, one can already anticipate that the zero-order reduced 
system will be inadequate for this model even with such simple coupling. Once the 
reference state x* is determined and z* computed, the zero-order reduced system is 
given, according to (|4.6|) . by 

_i + - ) {xi+i~Xi^2)- — + '^ -XyZ*. (5.5) 
a J a a'' 

This is equivalent to perturbing Fx by —a^XyJ*, which we expect to be small since 
X and y have zero mean in the uncoupled setting. In particular, we expect this 
perturbation to have only a small effect on the dynamics. However, in the multiscale 
dynamics it has been shown that a chaotic regime in the fast system can suppress 
chaos when coupled to the slow system [F', and this phenomenon is completely lost 
in the zero-order model. 

6. Numerical experiments 

Here we compare numerical results of the rescaled two-scale Lorenz '96 system 
with its corresponding reduced system. In particular we look at the ability of the 
reduced system captures some statistical quantities and how well it captures mean 
response to perturbations in the slow variables. 

In all parameter regimes considered, we have a slow system consisting of twenty 
variables (iVa; = 20) coupled with a fast system of eighty variables (iYy = 80). We use 
a fourth order Runge-Kutta method with timestep dt = e/10 in the multiscale system 
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and dt — 1/10 in the reduced system. To compute the mean response, an ensemble of 
10^ points is sampled from a single trajectory which has been allowed to settle onto 
the attractor. Using the translational symmetry of the Lorenz '96 system, we rotate 
the indices to generate an ensemble twenty times larger. 

On a modern laptop, the initial calculation to generate the reduced system for the 
Lorenz '96 system takes only a few minutes; once computed, numerical simulation of 
the reduced system is faster than the multiscale system by a factor on the order of e^^. 
Computing the mean response for a single forcing for 5 time units with a sufficiently 
large ensemble size (10^ trajectories) takes over an hour in the multiscale system with 
e = 10~^ but less than three minutes for the corresponding reduced system. 

6.1. Comparison of statistical properties of the two-scale and reduced 
systems 

In Section [3] we outlined the main requirements for correctly capturing the re- 
sponse of the two-scale system by its reduced model. Those were the approximation 
of joint distribution density functions (DDF) for slow variables, and the time auto- 
correlation functions of the time series. It is, of course, not computationally feasible 
to directly compare the 20-dimensional DDFs and time autocorrelations for all pos- 
sible test functions. However, it is possible to compare the one-dimensional marginal 
DDFs and simple time autocorrelations for individual slow variables, to have a rough 
estimate on how the statistical properties of the multiscale dynamics are reproduced 
by the reduced model. 

In figure 16.11 we compare the distribution density functions and autocorrelation 
functions of the slow variables. The DDFs are computed using bin-counting, and the 
autocorrelation function {xi {t)xi (t + s)), averaged over t, is normalized by the variance 
(xf). Results from three parameter regimes are presented, and in all three regimes the 
fast system is chaotic and weakly mixing {Fy = 12) and the coupling strength is chosen 
to be large enough (A^; = Aj, = 0.4) so that the multiscale dynamics are challenging to 
approximate. Of particular interest are timescale separations of £ = 10^^ and e = 10~^. 

First we consider a chaotic and strongly mixing slow regime {F^ = 16). Figures are 
presented for the timescale separation e — 10^^ only, because in this regime the picture 
is very similar for e = 10~^. We also consider a weakly chaotic and quasi-periodic slow 
regime {Fx = 8). In this regime, the coupled dynamics are more dependent on the 
timescale separation so we present results for both e — 10^^ and e — 10^^. Statistical 
quantities of other regimes, including regimes with more periodic behavior, have been 
presented in [4]. 

To more systematically compare DDFs for many parameter regimes, we introduce 
two metrics on the space of distributions. First is the Jensen-Shannon metric which is 
derived from the information-theoretic Kullback-Leibler divergence and is given 

by 

V2\J-oo \P{x) + q{x)J \p{x)+q{x)J J 

(6.1) 

where p and q are densities on distributions P and Q. The next metric we consider 
is the earth mover's distance |41| . which measures the minimum energy needed to 
move one DDF to another as though they were piles of dirt; the energy cost is the 
amount of 'dirt' times the distance it moved. One nice property is that the distance 
between two delta distributions 6a and 6p is simply ja — /?|. For distribution functions 
of one-dimensional random variables, the earth-mover's distance is the norm of the 



10 



The response of reduced models to small external perturbations 



F, = 16, £ = 0.1 
Distribution density Autocorrelation 



full multiscale system - 
zero-order closure 
first-order closure" 






full multiscale system — 




zero-order closure 




first-order closure — 











10 
Time 



i^^ = 8, e = G.l 



Distribution density 



Autocorrelation 



0.5 




full multiscale system — 
zero-order olosure 
first-order olosure — 


0.4 


1 \^ 


\ 


^0.3 

Q 




'A 


Q- 

0.2 


/ ^ 

I 




0.1 


/; 
i 


V 





3 -2 -1 
X 


1 2 : 




Distribution density 


0.5 




full multlsoale system — 
zero-order olosure 
first-order olosure — 


0.4 


/ V 






ft y 


\ \ 
A 1 
1 








0.2 






0.1 


'1 

■ '/ 
// 


\\ 

\s 
V 

V 





// 
y 




3-2-1 
X 


1 2 : 





full multiscale system — 




zero-order closure 




first-order closure — " 




/ ■^V 


/ / \^ 


\ / \ / 


t / / \\ 


\ / V / 


V • / / \ \ 











Autocorrelation 




Fig. 6.1. Distribution density and autocorrelation functions of slow variables. 



difference of the cumulative distributions 
mEM{P,Q) 



p{s)-q{s)ds 



dx. 



(6.2) 



Figure 16.21 shows distances between reduced systems DDFs and the corresponding 
multiscale slow variable DDFs. A variety of regimes is considered, with coupling 
parameters Xx,^y S [0.1,1], forcing parameters € {6,7,8,10,16} and Fy G {8,12,16}, 
and timescale separations e e {10~^, 10~^}. The data points are plotted with respect 
to coupling parameter A^. For each regime considered, the corresponding distances are 
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Fig. 6.2. Distances between DDFs of reduced and multiscale systems 



shown for both zero-order and first-order reduced systems. As the coupHng strength 
between fast and slow systems increases, it is apparently more difficult for the reduced 
systems to capture the correct slow dynamics of the multiscale system, as suggested 
by the linear best-fit curves. It should be noted that this correlation is slightly weaker 
when plotted against Xy, the coupling parameter for the fast system. However, the 
distribution densities of the first-order reduced system are consistently closer in both 
metrics than the zero-order system to the multiscale system. 

6.2. Mean state response to small perturbations 

In this section we examine the response of the mean state (x) of the slow variables 
(that is, h{x)=x in p.ip ) in the Lorenz '96 system to two simple types of small 
external forcing: 

1. Heaviside step forcing 



2. Ramp forcing 



ramp(^) 



if<<0, 
fH ift>0. 



ift<0, 

iframp ift>0. 



where djj and t^ramp are constant vectors. To compute the response of the mean state 
(a;), we generate an initial ensemble sampled from a trajectory that has been given 
sufficient time to settle onto the attractor. For each ensemble member we let evolve a 
short trajectory under the unperturbed dynamics as well as the under the perturbed 
dynamics, then we take the difference between these two trajectories and average over 
the entire ensemble. Here we consider forcing of the form v = aej , where a is constant 
and ej a standard basis vector in R^"^ . In the translation-invariant Lorenz '96 system, 
without loss of generality we need only consider a single such vector, say cq. Figure 
16.31 shows the mean response of the slow variables to small Heaviside forcing in the 
two-time rescaled Lorenz '96 system, where the magnitude of the forcing |t;H| is 1% 
of I /I averaged over the invariant distribution. 

Since the external forcing is sufficiently small to consider the response of the 
mean state approximately linear, we use the ideal response operator of Gritsun and 
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0.04 . - 




Fig. 6.3. Mean response of x to step forcing at node xn. Nx = '20, Ny = 80, Fx = 16, Fy = 
12, A^ = 0.4, Xy=OA, £ = 0.1 



Dymnikov (also see [TH51[51- [m[?7] ) by generating mean responses for several per- 
turbations and computing the linear best least-squares fit. This is a time-dependent 
matrix, and due to the symmetry of the Lorenz '96 system the dimensionality is 
reduced by one so that we have simply a time-dependent vector. With the ideal re- 
sponse operator, the response to a multitude of forcings can be readily estimated. We 
verify the nonlinearity in the actual response in Figure 16.41 which shows the growth 
in time of the relative error between the ideal response and the actual response to 
small Heaviside step forcing, averaged over several different forcings. We limit the 
plot to 5 time units because the large features of the Heaviside response are mostly 
fully developed by that time. Note in particular that the response is more linear in 
the reduced system, which is probably due to the fact that the Lyapunov exponents 
in the reduced system are much smaller than those in the two-scale system. 



Mean relative error of ideal response 




Multiscale system Reduced system 



Fig. 6.4. Nonlinearity of response: relative error ideal response vs actual nonlinear response. 




Fig. 6.5. Snapshots of the response operators for the response time T = 2 (left), and T = 5 
(right), Heaviside forcing. 



We now compare the ideal responses of the full and reduced systems. The snap- 
shots of the ideal responses for the two-scale and reduced models (as well as the linear 
response approximation, described in the next section) at times T — 2 and T — 5 are 
shown in Figures 16.51 and 16.61 for the Heaviside and ramp forcing, respectively. The 
response is captured accurately at the node which is directly forced (here Xn), but 
capturing the off-diagonal response as the perturbation propagates through the system 
is more difficult. Indeed, it seems that in the zero-order reduced systems the prop- 
agation speed is slightly faster than in the two-scale systems, but for the first-order 
reduced systems the response is well captured. 

It is more clear to see the quantitative differences in Figures 16.71 and 16.81 which 
show the relative distance between the responses as well as their cosine similarity 
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Response to ramp forcing. F^ — IG, e = 0.1 




5 10 15 20 5 10 15 20 



Response to ramp forcing. Fa; = 8, e = 0.1 




5 10 15 20 5 10 15 20 



Response to ramp forcing. Fx = 8, e = 0.01 




5 to 15 20 5 10 15 20 



Fig. 6.6. Snapshots of the response operators for the response time T = 2 (left), and T = 5 
(right), ramp forcing. 



II versus time for Heaviside step forcing and ramp forcing, respectively. Also 

shown in these figures are the linear responses computed using the reduced system 
statistics, as described in the next section. 

We observe that the first-order reduced system ideal response is a much closer 
approximation to the multiscale ideal response than the corresponding zero-order ideal 
response. In these regimes the relative error of the first-order response is limited to 
about 20% for the Heaviside forcing and less for the ramp forcing, while in the zero- 
order system the error is around 40% for the step forcing and 30% for ramp forcing 
at time t = 5. Remark that in the third plot for ramp forcing response in Figure [^751 
there is a small bump in the relative error shortly after the onset of forcing. This 
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Response to Heaviside forcing. F^ — W, e = 0.1 
Relative error Cosine similarity 



first- order ideal response 

zero- order ideal response 

first-order FDT response - ■ - 
zero -order FDT response — 



first- order FDT re 



Response to Heaviside forcing. = 8, £ = 0.1 
Relative error Cosine similarity 



first- order ideal response — 



;t-order FDT response ~ 
o-order FDT responaa - - 



Response to Heaviside forcing. Fx — 8, e = 0.01 
Relative error Cosine similarity 




Fig. 6.7. Comparing multiscale ideal response with reduced system ideal & quasi- Gaussian 
response operators for Heaviside forcing 



plot corresponds to a weakly chaotic regime {F^ — 8, Fy~ 12) with a large timescale 
separation (e = 10~^) in the multiscale system. In this regime the small nonlinear 
fluctuations of the multiscale system are relatively large compared to the ramp forcing 
for t near zero, so the relative error of the reduced system responses is large. 

6.3. Predicting the response of the two-scale system via linear response 
approximation of the reduced system 

Above in Section 16.21 we discussed the actual responses of the statistical mean 
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Response to ramp forcing. Fx = 16, e = 0.1 
Relative error Cosine similarity 



first- order ideal response 

zero- order ideal response 

first-order FDT response - ■ - 
zero -order FDT response — 



first- order FDT re 



Ramp forcing. F^ — S, e~0.1 
Relative error Cosine similarity 



first- order ideal response — 



?r FDT response" 
=r FDT responaa - - 




Response to ramp forcing. Fx = 8, £ = 0.01 
Relative error Cosine similarity 



first-order FDT re 



Fig. 6.8. Comparing multiscale ideal response with reduced system ideal & quasi- Gaussian 
response operators for ramp forcing 



states of both the two-scale and reduced models to small Heaviside and ramp forcings. 
For completeness of the study, we also attempt to predict the response of the mean 
state of the two-scale system via the quasi-Gaussian linear response approximation [I]- 
[3ll8l[9l[9l-[Tl ] [27 ] of the reduced system. In the quasi-Gaussian response approximation, 
the terms b{x) and Ba^x) in (|3.15p are replaced with the Gaussian approximations 
with same mean state and covariance matrices as in the actual dynamics. This, and 
the fact that h(x) = x in (|3.ip yields the following formula for the linear response 
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approximation of the mean state response: 




(6.3a) 



S{x)A{t)^ f RA{t-s)5f{s)ds, 

■JQ 

1 r 

i?A(i)= lim - / XAiT + t)){xAiT)~XAfdT^^\ 




(6.3b) 



where x and S are the mean state and covariance matrix of the corresponding unper- 
turbed systems (two-scale and reduced), computed as 



For large multiscale problems the mean response may be difficult to compute directly, 
since the large ensemble size needed for an accurate average is compounded by an al- 
ready large number of variables and small timestep discretization. In the case where 
the mean response of the slow variables is desired, one might prefer to compute the 
FDT response ()6.3|) using a time series from the reduced system for a "quick and 
dirty" approximation to the ideal response operator for the multiscale slow system. 
We show the accuracy of this FDT response approximation for the Lorenz '96 system, 
using a quasi-Gaussian approximation with time series data from the zero- and first- 
order reduced systems. The quasi-Gaussian response snapshots for the response times 
T = 2 and T = 5 are shown in Figures 16.51 and 16.61 for the Heaviside and ramp forcing, 
respectively. Qualitatively, the quasi-Gaussian response does capture the large fea- 
tures of the actual response, although most noticeable in these snapshots is the large 
exaggeration of the quasi-Gaussian response calculated from the first-order reduced 
system, which predicts a much larger off-diagonal response than what is observed. 
The possible reason for that is that the distribution densities of the first-order re- 
duced model are more strongly non-Gaussian than those of the zero-order reduced 
model, while the time autocorrelation functions are more weakly decaying (see Figure 
16. ip . It was observed previously in [27] that in these conditions the quasi-Gaussian 
linear response approximation tends to overshoot the off-diagonal response by a large 
margin. In other words, the better precision of the quasi-Gaussian linear response of 
the zero-order model is the result of mutual cancellation of the two errors: first one 
is the error in the distribution density of the zero-order reduced model (significantly 
more Gaussian than in in the multiscale dynamics) , while the second one is the error 
in the quasi-Gaussian linear response due to non-Gaussianity of the statistical state 
(less in the case of the zero-order model). 

The relative error and cosine similarity are measured against the multiscale ideal 
response and can be seen for Heaviside step forcing in Figure 16.71 and for ramp forc- 
ing in Figure 16.81 The ideal response of the first-order reduced system is clearly 




(6.4a) 




(6.4b) 
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the best of the four responses at capturing response in the slow variables. It is in- 
teresting, but perhaps not too surprising, that the least accurate estimate is given 
by the FDT response of the first-order reduced system. This should be expected 
since the quasi-Gaussian approximation is only valid for well-mixing systems whose 
distribution densities are close to Gaussian, which in particular is the case for the 
uncoupled Lorenz '96 systems in a chaotic regime F>8. However, such a system 
exhibits suppressed chaos when coupled to another chaotic systems, and the resulting 
distribution density will be far from Gaussian [5] . Since the first-order reduced system 
matches more closely the multiscale system, and the zero-order system will behave as 
an uncoupled Lorenz '96 system, the first-order system will be less chaotic and will 
be a poor candidate for the quasi-Gaussian FDT response. In fact, for non-chaotic 
regimes, as in the case of Fx = l, Fy = 12 where spatially periodic solutions emerge 
in the two-scale and first-order reduced systems, the long-time covariance matrix E 
will be singular, so the quasi-Gaussian response as presented will not be applicable. 
For further reading on blended FDT responses which might be more effective in these 
cases, see [8]. 

7. Summary 

In this work we studied the response to small external perturbations of multiscale 
dynamics and their reduced models for slow variables only. We elucidated a set of 
criteria for statistical properties of the multiscale and reduced systems which facili- 
tated similarity of responses of both systems to small external perturbations. It was 
shown that the similarity of marginal distribution densities of slow variables and their 
time autocorrelation functions controlled the similarity of responses to small external 
perturbations of both systems. 

Like in [1 , here we demonstrated that including a first-order correction term 
to a standard closure approximation for a nonlinear chaotic two-time system offered 
distinct improvements over the zero-order closure in capturing large-scale features 
of the slow dynamics. In particular, this reduced system was able to accurately 
capture the distribution density of solutions as well as the mean state response of the 
system to simple forcing perturbations. This correction term was relatively easy to 
generate, requiring only simple statistical calculations of the uncoupled fast system 
for an appropriate set of fixed parameters, and the resulting reduced system required 
much less computational resources than the corresponding multiscale system. 

Focusing on the mean state linear response of the slow variables, we showed that 
forcing perturbations in the reduced systems have similar responses as in the two-time 
system. Furthermore, we showed that using the unperturbed dynamics of the reduced 
systems for linear response prediction is also possible. However, in the parameter 
regimes we present here the first-order reduced systems are not rapidly mixing and do 
not follow a Gaussian distribution, but the zero-order reduced systems do have these 
properties, so this fluctuation-dissipation response is effective only using the zero-order 
system. A linear response method which takes into account the non-Gaussianity of 
the invariant statistical state (such as the blended response algorithm |8l-[TT|. based 
on the tangent map linear response approximation) is apparently needed to capture 
the response for strongly non-Gaussian dynamical regimes in reduced models. 

Here the linear response closure derivation and numerical results have been pre- 
sented only for the special case of linear coupling between slow and fast systems, 
but this derivation has been extended to systems with nonlinear and multiplicative 
coupling [5]. In future work we hope to extend similar results to these more general 
systems and to test the robustness of this method for application to a large variety of 
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problems. 
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